## Reproduction Script for:
## Mallinson, Daniel J. and Darrell Lovell
## "Race to the Top and the Diffusion of School 
## Turnaround Policy in the American States"
## Politics & Policy

rm(list=ls())

install.packages(c("dplyr", "cspp", "survival", "zoo", "mfx", "broom", "interplot", "gridExtra"))

library(foreign)
library(dplyr)
library(zoo)
library(lmtest)
library(mfx)
library(interplot)
library(gridExtra)
library(psych)
library(broom)

turnaround <- read.csv("turnaround_policy.csv")

year <- seq(min(turnaround$turnaround, na.rm=TRUE), max(turnaround$turnaround, na.rm=TRUE),1)

data <- as.data.frame(matrix(NA, nrow=1, ncol=7))
colnames(data) <- c("year", "state", "st.abb", "stateno", "state_fips", "state_icpsr", "adopt")

for(i in 1:nrow(turnaround)){
  a <- turnaround[i,]
  if(is.na(a$turnaround)){
    year <- seq(min(turnaround$turnaround, na.rm=TRUE), max(turnaround$turnaround, na.rm=TRUE),1)
    state <- a$statenam
    df <- as.data.frame(year)
    df$state <- a$state
    df$st.abb <- a$st.abb
    df$stateno <- a$stateno
    df$state_fips <- a$state_fips
    df$state_icpsr <- a$state_icpsr
    df$adopt <- 0
    data <- rbind(data, df)
  }else{
    year <- seq(min(turnaround$turnaround, na.rm=TRUE), a$turnaround,1)
    state <- a$statenam
    df <- as.data.frame(year)
    df$state <- a$state
    df$st.abb <- a$st.abb
    df$stateno <- a$stateno
    df$state_fips <- a$state_fips
    df$state_icpsr <- a$state_icpsr
    df$adopt <- 0
    df$adopt[df$year==a$turnaround] <- 1
    data <- rbind(data, df)
  }
}

data <- data[-1,]


## Add independent variables

# Citizen and government ideology
ideology <- read.csv("stateideology_v2018.csv")
data <- merge(x=data, y=ideology, by=c("state", "stateno", "year"), all.x=TRUE)

for(i in 1:50){
  q <- data[which(data$stateno==i),]
  state <- unique(q$state)
  q <- q[order(q$year),]
  q$citi6016 <- na.approx(q$citi6016, na.rm=FALSE, rule=2) #copies back and forward for trailing and leading NAs
  q$inst6017_nom <- na.approx(q$inst6017_nom, na.rm=FALSE, rule=2) #copies back and forward for trailing and leading NAs
  if(i==1){
    ideo.all <- q
  }else{
    ideo.all <- rbind(ideo.all,q)
  }
}

data <- ideo.all


#Neighbor Adoptions
neighbor <- read.csv("neighbors.csv")

df <- as.data.frame(matrix(NA, nrow=1, ncol=ncol(data)+2)) #Create blank dataset for appending
names(df) <- c(colnames(data), "neighbor_sum", "neighbor_prop")

vars <- c("state", "year","stateno","adopt")
adopts <- data[which(data$adopt==1),]
adopts <- adopts[vars]
for(j in 1:nrow(data)){
  use.obs <- data[j,]
  if(use.obs$year==min(data$year)){
    use.obs$neighbor_sum <- 0
    use.obs$neighbor_prop <- 0
    df <- rbind(df, use.obs)}
  else{
    obs.neighbor <- neighbor[which(neighbor$stateno == use.obs$stateno),]
    year <- as.numeric(use.obs$year)
    use.adopts <- adopts[which(adopts$year < year),]
    obs.neighbor$count <- 0
    for(k in 1:nrow(use.adopts)){
      obs.neighbor$count[use.adopts[k,3]==obs.neighbor$neighbor_no] <- 1
    }
    use.obs$neighbor_sum <- sum(obs.neighbor$count)
    use.obs$neighbor_prop <- mean(obs.neighbor$count)
    df <- rbind(df, use.obs)}
}

data <- df[-1,]

data$neighbor_per <- data$neighbor_prop*100 #Convert proportion to percentage

#Time Counter
data$time <- data$year - min(data$year)

#Cruz-Aceves and Mallinson (2020) Relative Ideology
df <- as.data.frame(matrix(NA, nrow=1, ncol=ncol(data)+1))
names(df) <- c(colnames(data),"relideo_cam")

vars <- c("stateno", "year","citi6016", "adopt")
adopts <- data[which(data$adopt==1),]
adopts <- adopts[vars]
for(j in 1:nrow(data)){
  use.obs <- data[j,]
  if(use.obs$year==min(data$year)){
    ideo.first <- data[which(data$year==min(data$year)),]
    ideo <- mean(ideo.first$citi6016)
    use.obs$relideo_cam <- 0
    df <- rbind(df, use.obs)}
  else{
    ideo.adopts <- adopts[which(adopts$year<use.obs$year),]
    old.adopts <- ideo.adopts[which(ideo.adopts$year!=max(ideo.adopts$year)),]
    new.adopts <- ideo.adopts[which(ideo.adopts$year==max(ideo.adopts$year)),]
    if(nrow(old.adopts)==0){
      adopters.ideo <- mean(new.adopts$citi6016)
    }else{
      adopters.ideo <- (mean(old.adopts$citi6016) + mean(new.adopts$citi6016))/2
    }
    use.obs$relideo_cam <- abs(use.obs$citi6016 - adopters.ideo)
    df <- rbind(df, use.obs)}
}

head(df)
df <- df[-1,] 

head(df)

data <- df

#Legislative Professionalism
squire <- read.csv("legislativeprofessionalism.csv")
names(squire)[1] <- "state"
names(squire)[3] <- "stateno"

data <- merge(data, squire, by=c("state", "year", "stateno"), all.x=TRUE, all.y=TRUE)

for(i in 1:50){
  r <- data[which(data$stateno==i),]
  state <- unique(r$state)
  r <- r[order(r$year),]
  r$legprof_squire <- na.approx(r$legprof_squire, na.rm=FALSE, rule=2) #copies back and forward for trailing and leading NAs
  if(i==1){
    squire.all <- r
  }else{
    squire.all <- rbind(squire.all,r)
  }
}

data <- na.omit(squire.all)

data$legprof_squire <- data$legprof_squire*100

#Race to the Top Dummy

data$racetop <- 0
data$racetop[data$year>2008] <- 1

#Demographics
demo <- read.csv("income_population_bea.csv")
names(demo)[1] <- "st.abb"

data <- merge(data, demo, by=c("st.abb", "state_fips", "year"),all.x=TRUE)

for(i in 1:50){
  s <- data[which(data$stateno==i),]
  state <- unique(s$state)
  s <- s[order(s$year),]
  s$inc <- na.approx(s$inc, na.rm=FALSE, rule=2) #copies back and forward for trailing and leading NAs
  s$population <- na.approx(s$population, na.rm=FALSE, rule=2)
  if(i==1){
    demo.all <- s
  }else{
    demo.all <- rbind(demo.all,s)
  }
}

demo.all$perinc <- demo.all$inc/demo.all$population

data <- demo.all

#Add log time
data$logtime <- log(data$time + .01)

#Add per pupil expenditures
pupil89 <- read.csv("stfis891b.csv")
pupil89 <- pupil89[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]
pupil89$SURVYEAR <- 1989

pupil90 <- read.csv("stfis901b.csv")
pupil90 <- pupil90[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]
pupil90$SURVYEAR <- 1990

pupil91 <- read.csv("stfis911b.csv")
pupil91 <- pupil91[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]
pupil91$SURVYEAR <- 1991

pupil92 <- read.csv("stfis921b.csv")
pupil92 <- pupil92[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]
pupil92$SURVYEAR <- 1992

pupil93 <- read.csv("stfis931b.csv")
pupil93 <- pupil93[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]
pupil93$SURVYEAR <- 1993

pupil94 <- read.csv("stfis941b.csv")
pupil94 <- pupil94[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]
pupil94$SURVYEAR <- 1994

pupil95 <- read.csv("stfis951b.csv")
pupil95 <- pupil95[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]
pupil95$SURVYEAR <- 1995

pupil96 <- read.csv("stfis961b.csv")
pupil96 <- pupil95[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]
pupil96$SURVYEAR <- 1996

pupil97 <- read.csv("stfis971b.csv")
pupil97 <- pupil97[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]
pupil97$SURVYEAR <- 1997

pupil98 <- read.csv("stfis981b.csv")
pupil98 <- pupil98[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]
pupil98$SURVYEAR <- 1998

pupil99 <- read.csv("stfis991b.csv")
pupil99 <- pupil99[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]
pupil99$SURVYEAR <- 1999

pupil00 <- read.csv("stfis001b.csv")
pupil00 <- pupil00[c("SURVYEAR", "FIPST", "STABR", "NCE13", "ADA")]
names(pupil00)[2] <- "FIPS"

pupil01 <- read.csv("stfis011b.csv")
pupil01 <- pupil01[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil02 <- read.csv("stfis021d.csv")
pupil02 <- pupil02[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil03 <- read.csv("stfis031b.csv")
pupil03 <- pupil03[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil04 <- read.csv("stfis041b.csv")
pupil04 <- pupil04[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil05 <- read.csv("stfis051b.csv")
pupil05 <- pupil05[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil06 <- read.csv("stfis061b.csv")
pupil06 <- pupil06[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil07 <- read.csv("stfis071b.csv")
pupil07 <- pupil07[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil08 <- read.csv("stfis081b.csv")
pupil08 <- pupil08[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil09 <- read.csv("stfis091b.csv")
pupil09 <- pupil09[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil10 <- read.csv("stfis101b.csv")
pupil10 <- pupil10[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil11 <- read.csv("stfis111a.csv")
pupil11 <- pupil11[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil12 <- read.csv("stfis121a.csv")
pupil12 <- pupil12[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil13 <- read.csv("stfis131a.csv")
pupil13 <- pupil13[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil14 <- read.csv("stfis141a.csv")
pupil14 <- pupil14[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil15 <- read.csv("stfis151a.csv")
pupil15 <- pupil15[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil16 <- read.csv("stfis161a.csv")
pupil16 <- pupil16[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil17 <- read.csv("stfis171a.csv")
pupil17 <- pupil17[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil18 <- read.csv("stfis181a.csv")
pupil18 <- pupil18[c("SURVYEAR", "FIPS", "STABR", "NCE13", "ADA")]

pupil19 <- pupil18
pupil19$SURVYEAR <- 2019

pupil <- rbind(pupil89, pupil90, pupil91, pupil92, pupil93, pupil94, pupil95,
               pupil96, pupil97, pupil98, pupil99, pupil00, pupil01, pupil02, pupil03,
               pupil04, pupil05, pupil06, pupil07, pupil08, pupil09, pupil10, pupil11,
               pupil12, pupil13, pupil14, pupil15, pupil16, pupil17, pupil18, pupil19)

names(pupil) <- c("year", "state_fips", "st.abb", "nce13", "ada")
pupil$nce13 <- as.numeric(pupil$nce13)
pupil$ada <- as.numeric(pupil$ada)
pupil$state_fips <- as.numeric(pupil$state_fips)

pupil <- pupil[which(pupil$state_fips <= 56),]

pupil$exp_pupil <- pupil$nce13/pupil$ada

data <- merge(data, pupil, by=c("year", "state_fips", "st.abb"))

#Exp Per Capita in $1,000s
data$exp_pupil10000 <- data$exp_pupil/10000

## Adjust pupil spending for inflation
year <- c(1989:2019)
cpi.values <- c(124, 130.7, 136.2, 140.3, 144.5, 148.2, 152.4, 156.6, 160.5, 163, 166.6,
                172.2, 177.1, 197.9, 184, 188.9, 195.3, 201.6, 207.3, 215.3, 214.5, 218.1, 224.9,
                229.6, 233, 236.7, 237.8, 240, 245.1, 251.1, 255.7)
cpi <- as.data.frame(cbind(year, cpi.values)) #https://www.minneapolisfed.org/about-us/monetary-policy/inflation-calculator/consumer-price-index-1913-

data <- merge(data, cpi, by="year")

data$exp_pupil10000_adj <- data$exp_pupil10000*255.7/data$cpi.values #Calculated in 2019 dollars

## Add Correlates of State Policy Variables
cspp <- read.csv("cspp_updated.csv")
cspp$st <- cspp$state <- cspp$state_fips <- cspp$state_icpsr <- NULL

for(i in 1:50){
  t <- cspp[which(cspp$stateno==i),]
  state <- unique(t$state)
  t <- t[order(t$year),]
  t$mathscore4th <- na.approx(t$mathscore4th, na.rm=FALSE, rule=2) #copies back and forward for trailing and leading NAs
  t$readscore4th <- na.approx(t$readscore4th, na.rm=FALSE, rule=2) #copies back and forward for trailing and leading NAs
  if(i==1){
    cspp.all <- t
  }else{
    cspp.all <- rbind(cspp.all,t)
  }
}

data <- merge(data, cspp.all, by=c("year", "stateno"))

#Create indicator of Republican state control
data$republican <- 0
data$republican[data$ranney4_control==0] <- 1

#Add Urbanization

urban <- read.csv("urban.csv")
data <- merge(data, urban, by=c("state", "year"), all.x=TRUE)

for(i in 1:50){
  u <- data[which(data$stateno==i),]
  state <- unique(u$state)
  u <- u[order(u$year),]
  u$pcturban <- na.approx(u$pcturban, na.rm=FALSE, rule=2) #copies back and forward for trailing and leading NAs
  if(i==1){
    data2 <- u
  }else{
    data2 <- rbind(data2,u)
  }
}

data <- data2

## Add indicator for RTTT Applicant
data$rttt_app <- 1
data$rttt_app[data$state=="Texas"] <- 0
data$rttt_app[data$state=="North Dakota"] <- 0
data$rttt_app[data$state=="Vermont"] <- 0

## Add indicator for RTTT winner
data$rttt_win <- 0
data$rttt_win[data$state=="Massachusetts"] <- 1
data$rttt_win[data$state=="New York"] <- 1
data$rttt_win[data$state=="Pennsylvania"] <- 1
data$rttt_win[data$state=="New Jersey"] <- 1
data$rttt_win[data$state=="Delaware"] <- 1
data$rttt_win[data$state=="Maryland"] <- 1
data$rttt_win[data$state=="Ohio"] <- 1
data$rttt_win[data$state=="Illinois"] <- 1
data$rttt_win[data$state=="Kentucky"] <- 1
data$rttt_win[data$state=="Tennessee"] <- 1
data$rttt_win[data$state=="North Carolina"] <- 1
data$rttt_win[data$state=="Georgia"] <- 1
data$rttt_win[data$state=="Florida"] <- 1
data$rttt_win[data$state=="Louisiana"] <- 1
data$rttt_win[data$state=="Colorado"] <- 1
data$rttt_win[data$state=="Arizona"] <- 1

## Export Analysis Dataset
write.csv(data, file="turnaround_analysis_data_rr.csv")

######################### Analysis ##########################################

data <- read.csv("turnaround_analysis_data_rr.csv")

data$RTTT <- 0
data$RTTT[data$year %in% c(2009:2011)] <- 1
data$postRTTT <- 0
data$postRTTT[data$year %in% c(2012:2019)] <- 1

################## Table 1 ##############################

describe(data[c("adopt", "racetop", "neighbor_prop", "relideo_cam", "republican", "pcturban",
                "legprof_squire", "exp_pupil10000_adj", "mathscore4th", "readscore4th")], skew=FALSE, na.rm=TRUE, range=TRUE)

################## Table 2 ##############################

model1a <- glm(adopt~ neighbor_per + relideo_cam + time, family=binomial(link="logit"), data=data)
model1b <- glm(adopt~ neighbor_per + relideo_cam + logtime, family=binomial(link="logit"), data=data)
lrtest(model1a, model1b)

## Diffusion Model
model1 <- logitmfx(adopt~ neighbor_prop + relideo_cam + racetop, data=data, clustervar1="state", atmean=FALSE)
model1
tidy(model1, conf.int=TRUE, conf.level=0.95)

#Odds Ratio
model1$oddsratio[,1]

#Lower CI
model1$oddsratio[,1] - 1.96*model1$oddsratio[,2]

#Upper CI
model1$oddsratio[,1] + 1.96*model1$oddsratio[,2]

#AIC
model1$fit

## NCLB Robustness Check
data$nclb <- 0
data$nclb[data$year %in% c(2002:2008)] <- 1
model1a <- logitmfx(adopt~ neighbor_prop + relideo_cam + nclb + RTTT + postRTTT, data=data, clustervar1="state", atmean=FALSE)
model1a
tidy(model1a, conf.int=TRUE, conf.level=0.95)

## RTTT Era Robustness Check
model1b <- logitmfx(adopt~ neighbor_prop + relideo_cam + turndate + afterturn, data=data, clustervar1="state", atmean=FALSE)
model1b
tidy(model1b, conf.int=TRUE, conf.level=0.95)


## Internal only model
model2 <- logitmfx(adopt~ republican + legprof_squire + exp_pupil10000_adj + mathscore4th + readscore4th + time, data=data, clustervar1="state", atmean=FALSE)
model2 
tidy(model2, conf.int=TRUE, conf.level=0.95)

#Odds Ratio
model2$oddsratio[,1]

#Lower CI
model2$oddsratio[,1] - 1.96*model2$oddsratio[,2]

#Upper CI
model2$oddsratio[,1] + 1.96*model2$oddsratio[,2]

#AIC
model2$fit

model4 <- logitmfx(adopt~ relideo_cam + racetop + legprof_squire  + exp_pupil10000_adj + racetop*exp_pupil10000_adj, data=data, clustervar1="state", atmean=FALSE)
model4
tidy(model4, conf.int=TRUE, conf.level=0.95)

model4$fit

#Robustness check with Cox model
model5 <- coxph(Surv(time, time2, adopt)~neighbor_prop + relideo_cam + ranney4_control + legprof_squire + exp_pupil1000_adj + mathscore4th + readscore4th + cluster(state), data=data,)
model5

model7 <- logitmfx(adopt~relideo_cam + legprof_squire + time + racetop*rttt_app, data=data, clustervar1="state", atmean=FALSE)
model7
tidy(model7, conf.int=TRUE, conf.level=0.95)
model7$fit

model8 <- logitmfx(adopt~relideo_cam + legprof_squire + time + racetop*rttt_win, data=data, clustervar1="state", atmean=FALSE)
model8

###################################### Figure 2 ##########################################

justadopts <- data[c("year", "adopt")]
justadopts <- aggregate(justadopts, by=list(justadopts$year), FUN=sum)
justadopts <- justadopts[c("Group.1", "adopt")]
names(justadopts)[1] <- "year"
justadopts$sum <- cumsum(justadopts$adopt)

png("figure2.png", height=6, width=7, res=600, units="in")
plot.new()
plot.window(xlim=c(1989,2020), ylim=c(0,50))
par(mar=c(2,4,0,0))
lines(justadopts$sum~justadopts$year)
axis(1, at=seq(1990,2020,5), labels=seq(1990,2020,5))
axis(2, at=seq(0,50,5), labels=seq(0,50,5), las=2)
title(ylab="Cumulative Count of Adopting States")
abline(v=2009)
text(2011, 40, "Race to the Top")
dev.off()

years <- seq(1989,2019,1)

for(i in 1:length(years)){
  if(i == 1){
    adopts <- data[which(data$adopt==1 & data$year==years[i]),]
    count <- sum(adopts$adopt)
    overtime <- as.data.frame(cbind(count, years[i]))
    names(overtime) <- c("count", "year")
  }else{
    adopts <- data[which(data$adopt==1 & data$year==years[i]),]
    if(nrow(adopts)==0){
      add <- cbind(count, years[i])
      names(add) <-  c("count", "year")
      overtime <- rbind(overtime)
    }else{
      adopts <- data[which(data$adopt==1 & data$year==years[i]),]
      count <- count + sum(adopts$adopt)
      add <- cbind(count, years[i])
      names(add) <-  c("count", "year")
      overtime <- rbind(overtime, add)
    }
  }
}

############################ Figure 3 ####################################

library(sandwich)
library(lmtest)
# https://www.r-bloggers.com/2021/05/clustered-standard-errors-with-r/

model9 <- glm(adopt~relideo_cam + legprof_squire + time + rttt_app*racetop, data=data, family=binomial(link="logit"))
summary(model9)

m9_coeffs <- coeftest(model9, vcov=vcovCL, cluster=~state)
m9_coeffs

#Code adapted from https://rpubs.com/milesdwilliams15/381372
ci <- .95
alpha <- 1-ci
z <- qnorm(1-alpha/2)
beta.hat <- coef(model9)
cov <- vcovCL(model9, cluster=~state)
z0 <- seq(min(model.frame(model9)[,"racetop"],na.rm=T),max(model.frame(model9)[,"racetop"],na.rm=T),1)
dy.dx <- beta.hat["rttt_app"] + beta.hat["rttt_app:racetop"]*z0
se.dy.dx <- sqrt(cov["rttt_app","rttt_app"] + z0^2*cov[nrow(cov),ncol(cov)] + 2*z0*cov["rttt_app",ncol(cov)])
upr <- dy.dx + z*se.dy.dx
lwr <- dy.dx - z*se.dy.dx

png("figure3.png", height=6, width=6, units="in", res=600)
plot.new()
par(mar=c(4,4,0,1.5))
plot.window(xlim=c(0,1), ylim=c(-1,15))
points(c(0,1),dy.dx, pch=20)
axis(1,c(0,1), c("Pre-RTTT","Post-RTTT"))
axis(2, seq(0,15,1), labels=seq(0,15,1), las=2)
segments(0,lwr[1],0,upr[1])
segments(1,lwr[2],1,upr[2])
abline(h=0, lty=2)
title(ylab="Effect of Being an RTTT Applicant")
dev.off()

###################################### Figure 4 ##########################################

model6 <- glm(adopt~ relideo_cam + racetop + legprof_squire  + exp_pupil10000_adj + racetop*exp_pupil10000_adj, data=data, family=binomial(link="logit"))
summary(model6)


a <- interplot(m=model6, var1="exp_pupil10000_adj", var2="racetop", hist=FALSE) +
  xlab("Race to the Top Era") +
  ylab("Marginal Effect of Net Expenditures per Pupil ($10,000s)") +
  theme_bw() +
  geom_hline(yintercept=0, linetype="solid")

b <- interplot(m=model6, var1="racetop", var2="exp_pupil10000_adj", hist=FALSE) +
  xlab("Net Expenditures per Pupil ($10,000s)") +
  ylab("Marginal Effect of Race to the Top") +
  theme_bw() +
  geom_hline(yintercept=0, linetype="solid")

png("figure4.png", height=5, width=10, units="in", res=600)
grid.arrange(a,b,nrow=1)
dev.off()

